%% 第七章：MIMO信道建模 - 模块化主仿真脚本
% MIMO信道建模与空间相关性分析
% 作者：周勇
% 日期：2024年

clc;
clear all;
close all;

%% 添加路径
addpath('../Common');

%% 获取颜色定义
colors = color_definitions();

%% 基本参数设置
carrier_freq = 2.4e9; % 载波频率 (Hz)
speed_kmh = 50; % 移动速度 (km/h)
speed_ms = speed_kmh * 1000 / 3600; % 转换为m/s
wavelength = 3e8 / carrier_freq; % 波长
max_doppler = speed_ms / wavelength; % 最大多普勒频移

% MIMO系统参数
num_tx_antennas = [2, 4, 8]; % 发送天线数量
num_rx_antennas = [2, 4, 8]; % 接收天线数量
antenna_spacing = 0.5; % 天线间距 (波长)

fprintf('=== MIMO信道建模 ===\n');
fprintf('载波频率: %.1f GHz\n', carrier_freq/1e9);
fprintf('移动速度: %d km/h\n', speed_kmh);
fprintf('最大多普勒频移: %.1f Hz\n', max_doppler);

%% 7.1 MIMO信道物理模型分析
fprintf('\n--- 7.1 MIMO信道物理模型分析 ---\n');

% 调用物理模型分析模块
results_physical = mimo_physical_models(...
    'carrier_freq', carrier_freq, ...
    'num_paths', 20, ...
    'angle_spread_tx', 30, ...
    'angle_spread_rx', 60, ...
    'mean_angle_tx', 0, ...
    'mean_angle_rx', 30, ...
    'antenna_spacing', antenna_spacing, ...
    'Nt', 4, ...
    'Nr', 4);

%% 7.2 MIMO信道容量分析
fprintf('\n--- 7.2 MIMO信道容量分析 ---\n');

% 调用MIMO信道分析模块
results_capacity = mimo_channel_analysis(...
    'SNR_dB', 0:5:25, ...
    'antenna_configs', num_tx_antennas, ...
    'num_realizations', 100, ...
    'carrier_freq', carrier_freq, ...
    'speed_kmh', speed_kmh);
    plot(singular_values.^2 / sum(singular_values.^2), 'o-', 'LineWidth', 2);
    grid on;
    xlabel('模式索引');
    ylabel('相对功率');
    title(sprintf('有效自由度: %.1f', effective_dof));
    
%% 7.3 天线阵列配置影响分析
fprintf('\n--- 7.3 天线阵列配置影响分析 ---\n');

% 调用天线阵列配置分析模块
results_array = antenna_array_config(...
    'SNR_dB', 0:5:25, ...
    'spacing_values', [0.1, 0.5, 1.0, 2.0], ...
    'array_configs', {'ULA', 'UPA', 'UCA'}, ...
    'carrier_freq', carrier_freq, ...
    'speed_kmh', speed_kmh, ...
    'Nt', 4, ...
    'Nr', 4);

%% 7.4 多普勒效应分析
fprintf('\n--- 7.4 多普勒效应分析 ---\n');

% 调用多普勒效应分析模块
results_doppler = doppler_effect_analysis(...
    'time_duration', 1, ...
    'sampling_rate', 1000, ...
    'speeds_kmh', [10, 50, 100], ...
    'carrier_freq', carrier_freq, ...
    'Nt', 4, ...
    'Nr', 4);

%% 结果保存与总结
fprintf('\n=== 结果保存 ===\n');

% 创建结果文件夹（在当前章节目录下）
if ~exist('results', 'dir')
    mkdir('results');
end

% 保存结果
save('results/mimo_modeling_results.mat', 'results_physical', 'results_capacity', 'results_array', 'results_doppler');

fprintf('结果已保存到 Chapter07_MIMO_Modeling/results/mimo_modeling_results.mat\n');

%% 总结
fprintf('\n=== 本章总结 ===\n');
fprintf('1. MIMO信道建模考虑了空间、时间和频率维度\n');
fprintf('2. 天线间距影响空间相关性，间距越大相关性越小\n');
fprintf('3. 阵列配置影响信道容量和实现复杂度\n');
fprintf('4. 多普勒效应导致信道时变，相干时间与速度成反比\n');
fprintf('5. 物理模型提供了真实的信道特性仿真\n');
fprintf('\n所有结果已保存到 Chapter07_MIMO_Modeling/results/ 文件夹\n');

fprintf('\n仿真完成！请查看生成的图形窗口。\n');

%% 辅助函数
function H = generate_physical_mimo_channel(Nt, Nr, num_paths, angle_spread_tx, angle_spread_rx, mean_angle_tx, mean_angle_rx, spacing)
    % 生成基于物理的MIMO信道
    H = zeros(Nr, Nt);
    
    % 角度转换为弧度
    mean_angle_tx_rad = deg2rad(mean_angle_tx);
    mean_angle_rx_rad = deg2rad(mean_angle_rx);
    angle_spread_tx_rad = deg2rad(angle_spread_tx);
    angle_spread_rx_rad = deg2rad(angle_spread_rx);
    
    for path = 1:num_paths
        % 随机角度
        angle_tx = mean_angle_tx_rad + angle_spread_tx_rad * randn();
        angle_rx = mean_angle_rx_rad + angle_spread_rx_rad * randn();
        
        % 随机复增益
        gain = sqrt(0.5) * (randn() + 1i * randn()) / sqrt(num_paths);
        
        % 发送端阵列响应
        a_tx = exp(-1i * 2 * pi * spacing * (0:Nt-1)' * sin(angle_tx));
        
        % 接收端阵列响应
        a_rx = exp(-1i * 2 * pi * spacing * (0:Nr-1)' * sin(angle_rx));
        
        % 累加路径贡献
        H = H + gain * a_rx * a_tx';
    end
end

function H = generate_correlated_mimo_channel(Nt, Nr, spacing, max_doppler)
    % 生成空间相关的MIMO信道
    % 发送端相关矩阵
    R_tx = zeros(Nt, Nt);
    for i = 1:Nt
        for j = 1:Nt
            R_tx(i,j) = exp(-0.5 * abs(i-j) * spacing);
        end
    end
    
    % 接收端相关矩阵
    R_rx = zeros(Nr, Nr);
    for i = 1:Nr
        for j = 1:Nr
            R_rx(i,j) = exp(-0.3 * abs(i-j) * spacing);
        end
    end
    
    % 相关信道矩阵
    R_sqrt_tx = sqrtm(R_tx);
    R_sqrt_rx = sqrtm(R_rx);
    H_uncorr = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
    H = R_sqrt_rx * H_uncorr * R_sqrt_tx';
end

function H = generate_array_mimo_channel(Nt, Nr, array_type, spacing)
    % 生成不同阵列配置的信道
    switch array_type
        case 'ULA'
            % 均匀线阵
            H = generate_correlated_mimo_channel(Nt, Nr, spacing, 0);
            
        case 'UPA'
            % 均匀面阵 (简化为两个ULA的组合)
            Nt_side = sqrt(Nt);
            Nr_side = sqrt(Nr);
            if Nt_side == floor(Nt_side) && Nr_side == floor(Nr_side)
                % 方形阵列
                H_horizontal = generate_correlated_mimo_channel(Nt_side, Nr_side, spacing, 0);
                H_vertical = generate_correlated_mimo_channel(Nt_side, Nr_side, spacing, 0);
                H = kron(H_vertical, H_horizontal);
            else
                H = generate_correlated_mimo_channel(Nt, Nr, spacing, 0);
            end
            
        case 'UCA'
            % 均匀圆阵
            H = generate_correlated_mimo_channel(Nt, Nr, spacing, 0);
            % 添加圆阵特有的相位关系
            radius = Nt * spacing / (2 * pi);
            for i = 1:Nr
                for j = 1:Nt
                    angle = 2 * pi * (j-1) / Nt;
                    phase_shift = exp(-1i * 2 * pi * radius * cos(angle));
                    H(i,j) = H(i,j) * phase_shift;
                end
            end
    end
end

function H_time = generate_time_varying_channel(Nt, Nr, time_axis, speed, carrier_freq)
    % 生成时变MIMO信道
    wavelength = 3e8 / carrier_freq;
    max_doppler = speed / wavelength;
    
    num_time_samples = length(time_axis);
    H_time = zeros(Nr, Nt, num_time_samples);
    
    % 初始信道
    H_time(:, :, 1) = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
    
    % 时变信道 (基于Jakes模型)
    for t = 2:num_time_samples
        time_diff = time_axis(t) - time_axis(t-1);
        
        % Jakes相关性
        correlation = besselj(0, 2 * pi * max_doppler * time_diff);
        
        % 生成相关的新信道
        innovation = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
        H_time(:, :, t) = correlation * H_time(:, :, t-1) + sqrt(1 - correlation^2) * innovation;
    end
end

function water_level = water_filling_mimo(eigenvalues, total_power, noise_power)
    % MIMO注水算法
    num_modes = length(eigenvalues);
    
    % 排序特征值
    [sorted_eig, sort_idx] = sort(eigenvalues, 'descend');
    
    % 寻找最优水位线
    water_level = 0;
    
    for k = 1:num_modes
        temp_level = (total_power + noise_power * sum(1./sorted_eig(1:k))) / k;
        temp_powers = max(0, temp_level - noise_power./sorted_eig(1:k));
        
        if all(temp_powers >= 0)
            water_level = temp_level;
        else
            break;
        end
    end
end